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Abstract Ionization fronts, the sharp radiation fronts 
behind which H/Hc ionizing photons from massive stars 
and galaxies propagate through space, were ubiquitous 
in the universe from its earliest times. The cosmic dark 
ages ended with the formation of the first primeval 
stars and galaxies a few hundred Myr after the Big 
Bang. Numerical simulations suggest that stars in this 
era were very massive, 25 - 500 solar masses, with 
H II regions of up to 30,000 light-years in diameter. 
We present three-dimensional radiation hydrodynami- 
cal calculations that reveal that the I-fronts of the first 
stars and galaxies were prone to violent instabilities, 
enhancing the escape of UV photons into the early in- 
tergalactic medium (IGM) and forming clumpy media 
in which supernovae later exploded. The enrichment of 
such clumps with metals by the first supernovae may 
have led to the prompt formation of a second gener- 
ation of low-mass stars, profoundly transforming the 
nature of the first protogalaxies. Cosmological radi- 
ation hydrodynamics is unique because ionizing pho- 
tons coupled strongly to both gas flows and primordial 
chemistry at early epochs, introducing a hierarchy of 
disparate characteristic timescales whose relative mag- 
nitudes can vary greatly throughout a given calculation. 
We describe the adaptive multistep integration scheme 
we have developed for the self-consistent transport of 
both cosmological and galactic ionization fronts. 
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1 First Light and the Rise of Primeval Galaxies 

State of the art numerical simulations reveal that the 
first stars in the universe formed in isolation in small 
cosmological halos (spheroidal bound clumps of dark 
matter and H/He gas) at redshifts z ~ 20 - 30, or 



200 - 400 Myr after the Big Bang (lAbel et all 120021: 
Bromm et al.ll2002l : lNa"kamura fc Umemural200ll) . Soon| 

thereafter, gravitational mergers congregated these star 
forming halos into the first primitive galaxies. Pre- 
liminary results from minihalo collapse models suggest 
that primordial (or Pop III) stars were very massive, 
25 - 500 M (jO'Shea fc Normanl l2007h . They cre- 
ated H II regions that were enormous in comparison 
to those in the Galaxy today, 2.5-5 kpc in radius 



(Whalen, Abel, fc Norman 2004; Ki tavama et al.ll2004 ; 



Alvarez et al.ll200a lAbel et al.ll2007l ). Primeval ioniza- 
tion fronts, or I-fronts, were key to the rise of the first 
protogalaxies for several reasons. First, Pop III stars 
lacked strong winds, so their H II regions determined 
their circumstellar environments. They regulated the 
flow of the first heavy elements from primordial super - 
novae (SNe) into the early IGM dGreif et al.ll2007il2010h 
and may even have resulted in the prompt formation of 
a second generation o f low-mass stars in SNe remnants 
( Whalen et al.ll2008bl ). Second, the minihalos forming 
the first stars themselves coalesced in small swarms 
due to cluster bias. Radiation fronts from one star in 
the cluster could engulf nearby halos, photoevaporating 
them to varying degrees and eithe r suppressing or al- 
lowing new star formation in them jShapiro et aL 2004 ; 
Iliev et al.|[2005l: ISusa fc Umemural 2006: Whalen et al 



2008al ISusa et al.l 120091: I Whalen et alj I20T0I) . Thus, 
cosmological I-fronts in part governed the number and 
character of stars that were swept up into the first 
primitive galaxies. 

Once primeval galaxies formed and sustained cycli- 
cal star formation, they too propagated I-fronts into 
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the IGM, beginning the process of cosmological reion- 
ization by which the cold, neutral, mostly featureless 
universe of the recombination era was gradually trans- 



and clusters of galaxies we observe todav ( Ricotti et al. 


2002. 


2008: 


Wise & Abel 


2008: 


Wise & Cenl2009). Re- 



cent numerical models have demonstrated that the I- 
fronts of both individual Pop III stars and protogalax- 
ies formed H2 in their outer layers, a key coolant at 
high rcdshifts that may have seeded the formation of 



additional structure in the early universe (jRicotti et al 
200lHWhalen fc Normanll2008bl) . High-redshift I-fronts 
thus influenced the rise of structure on a variety of spa- 
tial scales in the early cosmos. 

In the Galaxy today, I-fronts from massive stars in 
the molecular cores of giant molecular clouds (GMCs) 
are thought to trigger new star formation by perturb- 
ing nearby cores o n the verge of gravitation al collapse 



(jPale et alj l2007t iGritschneder et alj liooj) . In this 



manner, successive fronts can set off a chain of star for- 
mation throughout the cloud. On the other hand, dy- 
namical instabilities in the fronts may degenerate into 
turbulent flows that support gas against collapse, off- 
setting triggered star formation. Yet, these same in- 
stabilities can form clumps that may later collapse into 
stars. How ionized flows regulate star formation in the 
Milky Way today remains to be fully understood. 



2 Ionizing Radiation Hydrodynamics in 
ZEUS-MP 

We perfor m our ionization front calculations with 
ZEUS-MP jWhalen fc Normanlliooj) . a massively par- 
allel Eulerian astrophysical hydrodynamics code that 
solves the equations of ideal fluid dynamics^: 



dp 
dt 
dpvj 
dt 
de 

dt 



-V • (pv) (1) 
-V • {pViw) - Vp - pV<S> V Q (2) 
-V • (ev) - pV ■ v - Q : Vv. (3) 



Here, p, e, and the w 2 ; are the mass density, in- 
ternal energy density, and velocity at each mesh 
point and p = (7 — l)e and Q are the gas pres- 
sure and the von Neumann-Richtmeyer artificial vis- 
cosity tensor. ZEUS-MP evolves these equations 
with a second-order accurate monotonic advection 

1 http: //lea. ucsd.edu / portal /codes / zeusmp2 



scheme in one, two, or three dimensions on Cartesian 
(XYZ), cylindrical (ZRP), or spherical polar (RTP) 
coordinate meshes. Our augmented version of the 
publicly-available code self-consistently coupl es primor- 
dial gas chemistry ( Whalen fc Normanll2008a ) and mul- 
tifrequency photon-conserv ing UV radiative transfer 
( Whalen fc Normanl [2008b]) to fluid dynamics for the 
radiation hydroynamical transport of cosmological I- 
fronts. 

2.1 Primordial H and He Chemistry 

We evolve H, H+, He, Hc+, Hc 2 +, H , H+, H 2 , and 
e~ with nine additional continu ity equations and th e 
noncquilibrium rate equations of lAnninos et al.l (jl997j) : 



dpi 
dt 



-V • ( Pi v) 



EE 

j k 



Pjk(T)pjPk -^njpj, (4) 



where f3jk is the rate coefficient of the reaction between 
species j and k that creates (+) or destroys (-) species 
i, and the K/j 3X6 the radiative rate coefficients. We 
assume that the species share a common velocity dis- 
tribution. Mass and charge conservation, which are not 
guaranteed by either chemical or advective updates, are 
enforced each time the fluid equations are solved. 

Microphysical heating and cooling due to photoion- 
ization and gas chemistry is coupled to the gas energy 
density by an isochoric update that is operator-split 
from updates to the fluid equations: 



r-A, 



(5) 



where T is the cumulative heating rate due to photons 
of all frequencies and A is the sum of the cooling rates 
due to collisional ionization and excitation of H and He, 
recombinations of H and He, inverse Compton scatter- 
ing (IC) from the CMB, bremsstrahlung emission, and 
H2 cooling. 

2.2 Radiative Transfer 



Our photon-conserving UV transpo rt (jMellema et al 
2006UAbel. Norman, fc Madaulll999h . which is distinct 
from the flux-limited diffusion native to the public re- 
lease of ZEUS-MP, solves the static approximation to 
the equation of transfer in flux form to compute ra- 
diative rate coefficients for the reaction network at ev- 
ery point on the coordinate mesh. As currently im- 
plemented, our code can transport photons from a 
point source centered in a spherical grid or in plane 
waves along the x or z-axes of Cartesian or cylindri- 
cal meshes. To model Pop III spectra, we discretize 
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the blackbody photon emission rates of a star with 40 
uniform bins from 0.755 to 13.6 cV and 80 logarithmi- 
cally spaced bins from 13.6 eV to 90 eV, no rmalizing: 
them by the total ionizing photon rates by ISchaererl 
(2002). The radiative reactions in our models arc listed 
in Table 1 of IWhalen fe Normanl |2008bf ). We calcu- 



luate H2 photodissociation rates along rays parallel to 
the direction of radiation flow using sclf-shiclding func- 
tions modified for t herma l broadening as prescribed by 
Draine fe Bertoldil (119961 ) to approximate the effects of 
gas motion. They are show n in equations 9 and 10 of 
Whalen fe Normanl (|2008bl ). 



2.3 Radiation Forces 

Radiation pressure due to ionizations occurs at the I- 
front itself and in recomb ining gas in the H II region. 
Whalen fe Norman ( 20061 ) found that the acceleration 
of fluid elements at the front was large but momentary, 
and only alters the velocity of the front by 1 - 2 km 
s . Momentum deposition within the H II region is 
only prominent where gas is very dense, like at the cen- 
ter of a minihalo being evaporated by its star. There, 
rapid successive cycles of ionization and recombination 
can impart radiation forces to the gas that are hun- 
dreds of times the strength of gravity at early time s 
(lower left panel of Figure 1 in lKitavama et al.ll2004 ). 
As these forces propel gas near the center of the halo 
out into the H II region, its densities and recombina- 
tion rates fall, so more ionizing photons from the star 
reach the I-front. This higher flux results in I-fronts 
that are faster than when such forces are not included. 
In UV breakout from the first star-forming halos, radi- 
ation forces speed up the D-typc I-front by 10 - 20%. 
However, this effect is transient: after the internal re- 
arrangement of gas deep within the H II region dilutes 
its interior, radiation forces there sharply fall (see also 
Krumholz fc Matznerl 120091: iDrainei 12010; lLopez et af 
201(1 We describe the timescales on which momen- 
tum updates are performed below. 

2.4 Adaptive Subcycling 



Their relative magnitudes can seamlessly evolve during 
a single application. For example, when an I-front prop- 
agates through a medium, photoheating times are often 
smaller than Courant times, and chemical time scales 
are usually shorter than either one. On the other hand, 
fossil H II regions can cool faster than they recombinc, 
so cooling times become shorter than chemical times. 
The key to solving all three processes self-consistently is 
to evolve each on its own timescale without restricting 
the entire algorithm to the shortest of the times. To 
successfully deal with both I-fronts and relic H II re- 
gions, an algorithm must adaptively reshuffle the time 
scales on which the three processes are solved. Implicit 
schemes are sometimes applied to stiff sets of differen- 
tial equations like those in our model because they are 
unconditionally stable over the Courant time. However, 
accurate I-front transport in stratified media often re- 
quires restricting updates to both the gas energy and 
fluid equations to photoheating times in order to cap- 
ture the correct energy deposition into the gas, and 
linear system solves over such short time scales would 
be prohibitive in more than one or two dimensions. En- 
forcing photon conservation in implicit schemes can also 
be problematic. 

We instead subcycle chemical species and gas energy 
updates over the minimum of the chemical and heat- 
ing/cooling times until the larger of the two has been 
crossed, at which point we perform full hydrodynam- 
ical updates of gas densities, energies, and velocities. 
These times are global minima for the entire grid. The 
chemical times arc defined in terms of electron flow to 
accommodate all chemical processes rather than just 
ionizations or recombinations. Adopting the minimum 
of the two times for chemistry and gas energy updates 
enforces accuracy in the reaction network when t c h e m 
becomes greater than the, such as in relic H II regions. 
Our adap tive subcycling scheme is d escribed in greater 
detail in (I Whalen fc Normanl [2008ah . 



3 Ionization Front Instabilities 



A hierarchy of highly disparate characteristic time 
scales arises when gas dynamics, radiative transfer, and 
primordial chemistry are solved in a given application. 
The three governing times are the Courant time, the 
chemical time 



iche 



0.1 



0.001n ff 



and the photohcating/cooling time 



the 



0.1- 



^ht/cool 



(0) 



(7) 



When a point source of monoenergetic ionizing pho- 
tons appears in a uniform H and He medium, they are 
not free to simply stream into the interstellar medium 
(ISM) because of the large cross section to photoion- 
ization in the gas. Instead, they emerge as a sharp 
wall of radiation, behind which gas is almost fully ion- 
ized and heated to a fewxlO 4 K and beyond which it 
is completely neutral. The thickness of the I-front sep- 
arating the neutral and ionized gas is a few mean free 
paths of the ionizing photons through the neutral gas. 
The initial propagation of the I-front through the gas 
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monoenergetic: 



10 5 Kblackbody: 



quasar: 




secondary ionizations by photoelectrons 




Fig. 1 The structure of an ionization front. Xe denotes ion- 
ized fraction and T is gas temperature. Top: monoenergetic 
UV photons just above the ionization edge of H. Center: a 
10 5 K Pop III star blackbody spectrum. Note the signif- 
icant broadening of the front. Bottom: power-law quasar 
x-ray source spectrum. Note that in the first two cases that 
temperature tracks the ionization fraction, but not in the 
quasar case. Here, hard photons penetrate well ahead of 
the front, strongly heating but not significantly ionizing the 
gas. 



Fig. 2 Thin-shell instability in a primordial I-front medi- 
ated by H2 cooling. 



is supersonic in that the ionization wave advances far 
more quickly than any hydrodynamic response by the 
hot gas in its wake. This is known as an R-type I-front. 
As the front recedes from the star it slows, partly due 
to geometric dilution of the ionizing flux and partly be- 
cause recombinations in the ionized gas remove photons 
that would otherwise have gone on to advance the front. 
If the medium is static (no gas motion allowed), the I- 
front comes to a halt when the H II region it bounds en- 
closes enough ionized gas that the global recombination 
rate equals the emission rate of source photons. The 
radius at which this occurs is known as the Stromgren 
radius. 

In reality, as the I-front slows to the Stromgren ra- 
dius the pressure wave that has built up in the hot 
ionized gas behind it overtakes it and passes through 
it, steepening into a shock as it does so. Thereafter, 
the H II region expands subsonically with respect to 
the sound speed in its interior, which is still supersonic 
with respect to the ambient neutral gas. As it expands, 
it bypasses the static Stromgren radius limit because 
its interior densities and therefore recombination rates 
fall, allowing source photons to continue to stream to 
the front. As the H II region grows it pushes aside more 
and more ambient gas which accumulates on its surface, 
causing the shock to detach from and move ahead of the 
I-front. This is known as a D-typc I-front. 

3.1 Thin-Shell Instabilities 

If this shocked neutral shell can radiativcly cool, it 
will collapse into a cold dense layer that is prone to 



fragmentation v ia Vishniac, or thin-shell, instabilities 
( Vishniadll983f ). Unlike similar features in the dense 
shells of astrophysical blast waves, these erupt into very 
violent instabilities because ionizing radiation oppor- 
tunistically escapes thr ough the cracks in the shell in 
jets ahead of the front (IGarcia-Segura &: Franco! Il996 ; 



Whalen fc Normanl l2008al) . These jets rapidly crum- 
ple into complicated flows that form dense clumps that 
can persist for the life of the star. In the Galaxy today, 
metals and dust radiatively cool the shell and trigger 
such instabilities. Cooling by atomic H and He lines 
in the shocked gas of cosmological I-fronts by itself is 
too weak to incite such instabilities. However, they can 
instead form by H2 cooling. 

Primordial I-fronts are driven by hard UV sources 
whose photons have a variety of mean free paths 
through neutral gas. Collectively, they broaden the 
front, which leads to ionized fractions of 10% and tem- 
peratures of 2000 - 3000 K in its outer layers as we 
illustrate in Figure 1. These are ideal conditions for 
gas phase catalysis of H2 by the H~ and H^ channels: 



H 



H~ +7 



(8) 



H + H + -> HJ + 7 H+ + H -> H 2 + H+ (9) 

Although much less efficient than metals, H2 is a far 
more effective coolant than H and He lines at 3000 K; 
sandwiched between the front and the shocked neu- 
tral shell, it cools the base of the shell sufficiently to 



Fig. 3 The shadow instability. 



trigger thin-shell overstabilities that erupt into violent 
radiation-driven instabiliti es, as we show in Figure 2 
( Whalen fc Normanll2008b ). In this simulation, we cen- 
ter a 10 5 K blackbody UV source approximating a 120 
M Q Pop III star in a H/He gas profile similar to that 
of a cosmological minihalo. When the front becomes 
D-type, instabilities due to H 2 cooling can be seen to 
appear at the base of the dense shell. More models 
are required to determine if the turbulent gas clumps 
survive until the death of the star and are prone to 
gravitational collapse. 

3.2 Shadow Instabilities 

In Figure 3 we show the breakout of a second kind 
of radiation hydrodynamical instabi lity in primor- 
dial I-fronts, the shadow instability ( Williams! 1999 ; 
Whalen &: Normanll2008bl ). These form when density 
perturbations are advccted through an R-type front. 
The perturbation does not trap the front, it merely dim- 
ples it. The walls of the dimple are exposed to a smaller 
ionizing flux because they are oblique to it; as a con- 
sequence, the walls of the dimple transition to D-type 
before the tip does, causing the tip to advance and elon- 
gate the dimple into a jet that abruptly breaks down 
into a violent instability. Such phenomena would have 
been unavoidable in the primeval universe, on scales 
ranging from UV breakout from the first star-forming 
halos to 100 kpc protogalactic I-fronts. The role of H2 
formation, thin-shell instabilities, and shadow instabil- 
ities in structure formation in early cosmological flows 
remains to be properly investigated. 

3.3 AOI Instabilities 

We show in Figure 4 the radiation front of a 25 M Q 
Pop III star enveloping a nearby star-forming halo 




Fig. 4 Angle of incidence (AOI) instabilities in the ioniza- 
tion front of a Pop III star enveloping a neighbor halo. 



Spikes in the front can be seen 
all along the arc surrounding the halo, but are increas- 
ingly prominent with altitude above its axis. This is an 
examp le of an I-front instability described by Williams] 
( 20021 ) . in which photons that arrive at the front at 
any angle other than 90° cause it to be uncondition- 
ally unstable: an angle of incidence (AOI) instability. 
Detailed linear perturbation analysis predicts that the 
greater the angle of deviation of the photons from pure 
incidence, the faster the growth of the unstable mode, 
which true of the protrusions in Figure 4. In this exam- 
ple, the unstable modes have little effect on star forma- 
tion in the halo because those closest to its core have 
the smallest amplitudes and fail to puncture it. Dy- 
namical instabilities have also been found in I-fronts 
engulfing Galactic molecular cloud cores in numerical 
models, but gas in those si mulations could efficiently 
cool by fine-structure lines (jMizuta et al.l 120051 120061 ) 
so the unstable modes are likely thin-shell in origin. 
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